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In ferromagnetic Bose-Einstein condensates (BECs), the quadratic Zeeman effect controls mag¬ 
netic anisotropy, which affects magnetic domain pattern formation. While the longitudinal mag¬ 
netization is dominant (similar to the Ising model) for a negative quadratic Zeeman energy, the 
transverse magnetization is dominant (similar to the XY model) for a positive one. When the 
quadratic Zeeman energy is positive, the coarsening dynamics is driven by vortex-antivortex anni¬ 
hilation in the same way as the XY model. However, due to superfiuid flow of atoms, there exist 
several combinations of vortex-antivortex pairs in ferromagnetic BECs, which makes the coarsening 
dynamics more complicated than that of the XY model. We propose a revised domain growth 
law, which is based on the growth law of the two-dimensional XY model, for a two-dimensional 
ferromagnetic BEC with a positive quadratic Zeeman energy. 
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I. INTRODUCTION 

Domain growth and coarsening dynamics have been 
studied in a wide variety of systems [iHlll • When a sys¬ 
tem is quenched from a disordered phase to an ordered 
phase, the long-range order does not arise immediately. 
At first, locally-ordered small domains arise, which then 
grow with time to make global order through domain 
coarsening. Although there are different mechanisms to 
cause domain growth, in most cases, domain size I grows 
with time t as l{t) ~ , where iz is an scaling exponent. 

For example, iz = 1/3 in the two-dimensional (2D) con- 
seryed systems described by the Ising model (e.g., binary 
alloys and ferromagnets with uniaxial anisotropy) [iH^. 
If fluid flow contributes to domain growth, which is the 
case of binary fluids, the exponent changes depending on 
advection and viscosity. When the advective transport 
is negligible, diffusion dominates the coarsening dynam¬ 
ics. In that case, iz = 1/3 [^, which is the same as 
in the absence of flow. However, if the advective trans¬ 
port with little viscosity dominates over diffusion, the 
inertia of fluid becomes important in the coarsening dy¬ 
namics. In this case, domains grow faster than the dif¬ 
fusive case, and the exponent is iz = 2/3 When the 
system is described by vector fields (i.e., complex order 
parameters or multi-component order parameters), the 
dominant mechanism to cause domain growth is com¬ 
pletely different from those for the Ising model and bi¬ 
nary fluids. For example, the coarsening dynamics for 
2D vector fields is driven by the annihilation of vortex- 
antivortex pairs. The domain size, which is actually the 
characteristic length of the spatial structure of the field, 
grows as l{t) ~ for non-conserved n-component vec¬ 
tor fields in d-dimensional space, except for d = n = 2, 
namely, the 2D XY model. The domain growth law 
for the 2D XY model includes a logarithmic correction: 


l{t) ^ (t/lnt)^/^ [3-[ni- 

Magnetic domain patterns and their coarsening dy¬ 
namics are observed also in ferromagnetic Bose-Einstein 
condensates (BECs). Recent development in imaging 
techniques to observe magnetization profiles in ferro¬ 
magnetic BECs has enabled us to investigate the real¬ 
time dynamics of magnetization, such as spin texture 
formation, spin-domain coarsening, and nucleation of 
spin vortices (T2l - [l5l| . Those experiments have also 
motivated theoretical studies about configurations of 
Skyrmions and spin textures Em , magnetic domain 
formation [ll, [l^, and spin turbulence |20l - [2^ . Mag¬ 
netic anisotropy of a ferromagnetic BEC depends on the 
quadratic Zeeman energy, which can be controlled by ex¬ 
ternal fields. When the quadratic Zeeman energy is neg¬ 
ative, longitudinal magnetization is dominant, and thus 
the system is similar to the Ising model. In 2D ferromag¬ 
netic BECs with a negative quadratic Zeeman e nerg y or 
binary BECs, domain size grows as l{t) ~ [l9.1^. 

which has the same exponent iz = 2/3 as that for binary 
fluids in the inertial hydrodynamic regime. However, 
l{t) ~ in the absence of superfluid flow [l^. The 
difference in the exponents suggests that the superfluid 
flow has a strong influence on the coarsening dynamics. 

In this paper, we investigate the coarsening dynam¬ 
ics in 2D spin-1 ferromagnetic BECs with a positive 
quadratic Zeeman energy. When the quadratic Zee- 
man energy is positive, transverse magnetization is dom¬ 
inant, and the coarsening dynamics is caused by vortex- 
antivortex annihilation. The situation is similar to the 
2D XY model, however a crucial difference arises in the 
classihcation of vortices. Vortices in ferromagnetic BECs 
are classified by the winding number of spin current (di¬ 
rection of magnetization) and mass current (vorticity of 
superfluid flow). Thus, there are several combinations 
of vortex-antivortex pairs which cause pair annihilation 
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in ferromagnetic BECs. By contrast, in the XY model, 
there is only one combination of vortex-antivortex pairs. 
When several combinations of annihilation pairs exist, 
the coarsening dynamics is expected to be more compli¬ 
cated than that of the XY model. In other words, super¬ 
fluid flow has indirect effects on the coarsening dynamics 
through different combinations of vortex-antivortex anni¬ 
hilation. We will demonstrate the coarsening dynamics in 
ferromagnetic BECs by numerical simulations, and pro¬ 
pose a revised domain growth law, based on the growth 
law for the XY model. 

The rest of the paper is organized as follows. The 
decay of the vortex density, which is caused by vortex- 
antivortex annihilation, in ferromagnetic BECs is dis¬ 
cussed in Sec. im Numerical simulations illustrated in 
Sec. mil clearly show that superfluid flow affects the 
coarsening dynamics and that a revised growth law is 
needed for the case where there are several combinations 
of vortex-antivortex pairs. The revised law is proposed in 
Sec. |TVl Discussions and conclusions are given in Sec. El 


II. DOMAIN GROWTH LAW 

A. Interaction of vortices 


We consider a spin-1 BEC confined in the x-y plane 
under a uniform magnetic field applied in the z direction. 
For simplicity, we neglect the confining potential in the 
X and y directions. The mean-field kinetic energy and 
Zeeman energy are given by 

= f dr ^2 ( 1 ) 

Eq= f dr (2) 

d m=-l 

where 'bm('f') is the condensate wave function for the 
atoms in the magnetic sublevel m, M is an atomic mass, 
and q is the quadratic Zeeman energy per atom. Here, 
we neglected the linear Zeeman term because the linear 
Zeeman effect merely induces the Larmor precession of 
atomic spins and can be eliminated in the rotating frame 
of reference. The quadratic Zeeman energy is tunable 
by means of a linearly polarized microwave field and can 
take both positive and negative values . 

The interatomic interaction energy is given by 

Eint = \ j [contotirf + ci\f{r)\'^] , (3) 

where the number density and the spin density (local 


magnetization) are given by 
1 

ntot{r) = Y ( 4 ) 

m— — l 

1 

fvi'f') = 'Y/ ^ m{'^){Ev)mn'^ nif") , (5) 


respectively. Here, v = x,y,z and, Fx,y,z are the spin- 
1 matrices. The interaction coefficients are given by 
Co = 47rIi^(2a2-l-ao)/(3M) and ci = 47r/i^(a2 —ao)/(3M), 
where as is the s-wave scattering lengths of two colliding 
atoms with total spin S channel. For the condensate to 
be stable, cq needs to be positive. On the other hand, 
the sign of ci determines the magnetism: the condensate 
is ferromagnetic (antiferromagnetic or polar) for ci < 0 
(ci >0). In this paper, we consider ferromagnetic BECs 
(ci < 0). 

When the quadratic Zeeman energy is weak compared 
with the ferromagnetic interaction energy, the conden¬ 
sate is fully magnetized (|/| = ritot)- Since the order 
parameter for a fully-magnetized state in the direction 
(cos a sin /3, sin a sin /3, cos /3) is given by [H, [13 



the population in the m = 0 component becomes max¬ 
imum at j3 = 7r/2, whereas those in the m = 1 and — 1 
components become maximum at /3 = 0 and tt, respec¬ 
tively. As seen from Eq. the quadratic Zeeman effect 
enhances the population in the to = 0 state for g > 0 and 
those in the to = ±1 states for q < Q. Hence, the mag¬ 
netization arises in the x-y plane [P = tt/2) for g > 0, 
and in the -hz or —z direction (/3 = 0 or tt) for g < 0. 
The former case corresponds to the XY model and the 
latter the Ising model of the conventional ferromagnet. 
Although the magnitude of the spontaneous magnetiza¬ 
tion becomes smaller when the quadratic Zeeman energy 
is positive and comparable to the ferromagnetic interac¬ 
tion, the magnetization direction is still confined in the 
x-y plane. Since we are interested in vortex-antivortex 
annihilation, we consider g > 0 below. 

We first consider a single vortex and write its wave 
function in the polar coordinate whose origin is the center 
of the vortex core: (r, ip). We take (p = and a = 
OaP) in Eq. ([H), where and ctq are integers. For a 
symmetric vortex, /3 is a function of r and independent 
oi ip. At a distance from the vortex core, /3 = 7r/2 as 
discussed in the above. The wave function outside of the 
core is approximately written as 



( 7 ) 


Here, and a a determine the directions of mass flow 
and spin flow around the vortex, respectively. The super¬ 
fluid velocity (mass flow) and the spin superfluid velocity 
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(spin flow) of z component are written for a homogeneous 

^tot 


Vmass — /^tot, 

m— — l 

( 8 ) 

h ^ 

'*^spin = 2Mi i^z)mn ~ /''^-totj 

m— — l 

(9) 




FIG. 1. Evaluation of the energy impair of two vortices with 
vorticities and separated by a distance D. 



respectively. Substituting Eq. © into Eqs. ([5]) and (O, 
we see that the directions of mass and spin flows depend 
on (70 and cr^, respectively, as i^mass = o-^{h/M)'S/(p and 
<pin = -^Ufil2M)V^. 

The combination of and Ca also determines the vor¬ 
tex core structure. Though we use /? = 7 r /2 in Eq. d?)), /3 
changes around the center of the vortex so as to remove 
the singularity of the order parameter. When cr^ = ctq, 
the m = 1 component is independent of ip and only this 
component remains at r = 0. In this case, (3 takes 0 at 
r = 0 , which means the magnetization at the center is in 
the +z direction for = <Ja- Similarly, for vortices with 
CT 0 = —<7a, magnetization is in the —z direction at the 
center. When tr^ = 0 and Oa 0, (/j-dependent compo¬ 
nents cannot vanish in a fully-magnetized state. Thus, 
magnetization vanishes at the center for cr^ = 0. On the 
other hand, when cr^ 0 and CTq = 0 , all three com¬ 
ponents should vanish at the center. In the following, 
we consider only the elementary vortices that are sta¬ 
ble against splitting, that is, = 0 , ±1 and a a = ± 1 - 
The vortex of cr^ = 0 has no mass flow around itself 
and its core is not magnetized. Such a vortex is called 
polar-core vortex (PCV). When = ±1, the vortex 
has mass flow around its core and its core is fully mag¬ 
netized as well as the outside. Such a vortex is called 
Mermin-Ho vortex (MHV). Considering the combina¬ 
tion of CT 0 and a on we notice that there are two kinds 
of PCVs [(cr 0 ,cra) = (0, ±1)] and four kinds of MHVs 
= (± 1 ,± 1 )]. 

A vortex-antivortex pair is defined so that they can 
be pair-annihilated. For the case of a single-component 
BEC, two vortices with winding numbers with the op¬ 
posite signs are a vortex-antivortex pair. In the present 
case, i.e., a multi-component BEC, when two vortices 
have winding numbers with the opposite signs in all com¬ 
ponents, they can be annihilated as a vortex-antivortex 
pair. For the vortex expressed by Eq. o, the m = I, 0, 
and —1 components have the winding numbers cr^ — ctq,, 
atf,, and + oa, respectively. Thus, its antivortex is 
obtained by changing the signs of both and a a. In 
other words, vortices with {a^,aa) and are 

a vortex-antivortex pair. 


we estimate the kinetic energy of a single vortex as 

m 

^2 pR |•2^T 

^He ^0 ^ 

= c'Ariiog(^^), (10) 


where C = 'nh^ntot/M, and R and Rc are the vortex size 
(radius) and the radius of the vortex core, respectively. 
Although R is equal to the system size for a single vor¬ 
tex, it is the distance beyond which the field around the 
vortex is shielded if there are other vortices. A/i depends 
on and a a ■ Since the portions of the number densities 


for m = 1 , 0 , —1 are t, respectively, 




( 11 ) 


where izi = — aa, vq = and iz-i = + (ja- 

The interaction energy between two vortices at a dis¬ 
tance of D is given by V{D) = Epair{D) — E 1 — E 2 , where 
Apair is the energy of two vortices separated by a distance 
D, El and E 2 are the energies of single vortices with vor¬ 
ticities and iym\ respectively. The pair energy Apair 
is approximately given by the sum of contributions from 
two regions (see Fig.[T|). In the region oi D < r < R, the 
contribution is evaluated for a single (composite) vortex 
with vorticity ■ In the region of r < U, the 

contribution is just the sum of Ei and i? 2 . Then, the 
interaction energy is approximated by 


V{D) = CM 2 log ( ^ 




^0 ^0 




( 12 ) 

(13) 


The derivative of V(D) gives the force between the vortex 
pair, 


dV{D) CM 2 
dD ~ D 


(14) 


In this paper, we consider only MHVs, which are useful 
to investigate the effect of superfluid flow. Using Eq. 0 , 


which is an attractive force between a vortex-antivortex 
pair. Note that pair annihilation occurs only between a 
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vortex pair in which both cr^ and ctq, have the opposite 
signs, although the force is attractive (A^ < 0) between 
vortices with the opposite signs of even if they have 
the same sign of Ua- 


B. Coarsening dynamics 


Here, we assumed a uniform number density: Vritot = 0. 
The kinetic energy in this formulation is written as 


-Ekin = 


AM 
, Mntot 


dr 


(V/.)2 + (V/,)2 + [VJ.f 


dr vi 


(19) 


We assume that the attractive force between a vortex- 
antivortex pair is balanced with a friction (resistive) force 
Tfric when the pair vortices move toward each other. The 
friction causes energy dissipation. When a vortex moves 
at speed u, the energy dissipation rate is written as 

dE 

7, ~ (1^) 

dt 


The dynamics of a spinor BEC is well described with 
the time-dependent multi-component Gross-Pitaevskii 
(GP) equation, and we phenomenologically introduce an 
energy dissipation into the GP equation [23, Ull : 




- d-it) + qm^ + contot{r,t) 




+ E mn ^ n 0; (16) 

n= —1 i^—x,y,z 


where P expresses energy dissipation. 

In order to discuss the energy dissipation that is caused 
by the friction force, we employ the hydrodynamic equa¬ 
tion, which is derived in the low-energy limit [ilili. In 
this limit, the BEG is fully magnetized, i.e., |/| = ntot 7 
and the physical quantities that describe the dynamics 
of ferromagnetic BEGs are the normalized spin vector 



“^tot 


(17) 


We divide the energy dissipation into two parts, which is 
written as 


dE 

dt 


dE„ 


dt 


dt 


dEmag dEf[ OW 

dt dt 




/ 5E dvx 
V<5z;a, dt 


5 fy 

SE dvy 

SVy dt 


Sf. dt j 


( 20 ) 

( 21 ) 

( 22 ) 


where t^mass = {vx,Vy). We assume that a vortex keeps 
its shape, i.e., the profiles of t^mass and / around its core, 
when it moves. The contribution to dEma,g/dt arises from 
the change in direction of local magnetization. Since the 
profiles of t^mass and / are conserved, the coupling be¬ 
tween Uinass and f gives no contribution to dE^^^g/dt. 
Neglecting the energy contributions from the vortex core, 
we only need to consider the hydrodynamic equation in 
the outside region of the vortex core. Then, Eq. (fT^ with 
fz = dfz/dt = Vfz = 0 leads to 


dfx 

dt 

dfy 

dt 


dt 


h 


1 -h r2 2M 
p h 
1 -h r2 2M 

V(V-t; 


2Mr 


fyif-^^f) 


vVy 


(23a) 

(23b) 

(23c) 


where the coupling terms with (t^mass ■ V)/ are dropped. 
Since we take fz = 0, fx = cos a and fy = sin a. From 
Eqs. (I23al) and (I23bl) . we have 


and the superfluid velocity t^mass- The equations of mo¬ 
tion for them are derived straightforwardly from the GP 
equation (HU) [ii,[ii,[2ii, and the resulting equations of 
motion are written as 


dt 


1 + T^ 

P 


/ X BeS ("^mass V)/ 

1 


1 -hP^ 


H 

/x 


:/ X Beff ('^mass V)/ 


fl2 

Bfiff = -rvyV^/ - qfzZ, 


2M 

TlkT ^ ^ 

-^■^'t^mass = - — 

dt 2ntotr 


V [V ■ 

ass)] 


+ /X 


dt 


(18a) 

(18b) 

(18c) 


da _ ^ dfy ^ df:z _ 

~ Jy ~ 


r 


-V^a. 


dt dt dt l + T'^2M 
Substituting Eqs. (I23al) and (I23bl) into Eq. (|2T|) gives 


(24) 


^ntot(l + P^) 


dr 


dfx 

dt 


dt 


hutotii+T"^) 


(25) 


where we used f^+fy = 1- Substituting Eq. dZD into 
Eq. ([5]), we have 


= —V((>, 
M 


(26) 









































5 


where ij) = a^(p. Equations. (I23c|) and (|26)) lead to 


dt 2MT ^ 


(27) 


Using Eqs. (I23c|) . (E51) . and (E71) . we rewrite Eq. (1^ as 


where to is an integration constant. Employing the vor¬ 
tex density p = 1 /^^ and the maximum vortex density 
Pc = l/7?cj '"^6 rewrite Eq. (IMll as 


^ 1 logjpclp) - 1 

AA p 


(35) 


— Mfltot I dv l^mass * 

= -2hntotT j dr . 


(28) 


Suppose that a vortex keeps its shape when it moves 
in the x direction at speed u: a{r) = f(x — ut,y) and 
'/'(’’) = 9(a ~ where / and g are functions ex¬ 

pressing their profiles. Then, (da/dt)^ = v?{da/dx)'^ 
and {d(j)/dt)‘^ = v? {d(j)/dxY. Similarly, for a vortex 
moving in the y direction, (da/dt)"^ = u^{da/dy)'^ and 
{d(j)/dt)‘^ = u‘^{d(j)/dy)‘^. The averages of them result in 


da 


dt 


d(j) 


^ =^(Va)^ =^(V</>)^. (29) 


dt 


Combining Eqs. (051) and (1^ and using Eq. (0^1) . we 
estimate the energy dissipation as 



= (30) 


Here, 

Mric = (1 + T^Wi - al, (31) 

where we used a = Ua^p. The friction force is estimated 
by comparing Eqs. m and O: 

2A(f / R \ 

Ffnc = -^CA/fric log ( j (32) 

In order to investigate the growth of characteristic do¬ 
main size we apply the discussion of the coarsening dy¬ 
namics in the 2D XY model P,[l,|^[ill, where we expect 
^ ~ i? ~ D and u = d^/dt. Equating the characteristic 
force between a vortex pair Epair oc A/ 2 /^ with the char¬ 
acteristic friction force Emc oc (A/fric/E) log(^/Ec)d^/dt, 
and rearranging terms, we have 

(1:) § ='*■ (=*=•> 


where A is a constant that depends on the dissipa¬ 
tion rate and the characteristics of vortices as A oc 
—EA^/A/fric- Integrating of Eq. (1^ gives 





= 2A{t-to), 


(34) 


The number of vortices in ferromagnetic BECs is ex¬ 
pected to yield Eq. (1551) . which is the same as the growth 
law for the XY model The difference between the 
XY model and ferromagnetic BECs is contained in factor 
A, which includes the information about vortices (Afric 
and A/ 2 ). Actually, the hydrodynamic equation with 
'I’mass = 0 corresponds to the XY model when ~ 0 
(namely, in a positive -9 case). If t^mass = 0, we just drop 
the Eflow terms in the above discussion, and then obtain 
the same equation as Eq. dSSl), although the factor A is 
different from that of the above case. 


III. NUMERICAL SIMULATIONS 

We perform numerical simulations by means of the dis¬ 
sipative GP equation (I16|) and the hydrodynamic equa¬ 
tion (ITSl) . The advantage of the hydrodynamic equa¬ 
tion is that the superfluid velocity Umass can be elimi¬ 
nated easily in simulations, which enables us to investi¬ 
gate what effects the superfluid flow has on the coarsen¬ 
ing dynamics. Note that MHVs introduced in Sec. Ill Al 
have both mass flow and spin flow around their cores. 
However, if we take Umass = 0 in simulations, the degrees 
of freedom of mass flow are eliminated. Then, the hydro- 
dynamic equation reduces to the equation of motion of 
magnetizations, and what we call MHV in the discussion 
below becomes merely a spin vortex around which only 
the spin current circulates. In such a case, the index 
is meaningless, and there exists only one combination of 
a vortex-antivortex pair; ctq, = I and —I. 

In the simulations, the mass of an atom is given by a 
typical value for a spin-I ®^Rb atom: M = 1.44 x 10“^^ 

kg^_The total number density is taken as ritot = 

whuPn^Y:, with 7130 = 2.3 x 10^^ cm“^ and d = I /rm. 
The quadratic Zeeman energy is set to be q/h = 10 
Hz. The dissipation rate is given by a typical value 
E = 0.03. Especially in the GP simulation, the sys¬ 
tem is in quasi-two dimensions: The wave function in 
the normal direction to the 2D plane is approximated 
by a Gaussian with width d. Interaction parameters are 
taken as consu/h =1.3 kHz and citisd = —59 Hz. The 
value of Cl that we take here is ten times larger than a 
typical value of a spin-I ®^Rb atom, which prevents the 
production of PGVs. 

Initial states are given by randomly located four kinds 
of MHVs. The number of vortices of each kind is equal. 
Open boundary conditions are imposed on 256 ^m x 256 
pm systems. The total number of vortices at first are 256, 
which implies that the average distance between vortices 
is about 16 ^m. Snapshots of the transverse and longitu¬ 
dinal magnetizations, the vorticity of mass flow, and the 
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FIG. 2. (Color online) Snapshots of the transverse 
(a,Tcta,n{fy/fx)) and longitudinal (fz) magnetizations, the 
vorticity of mass flow (V x Vmass), and the positions of vortex 
cores at time t = 1 s (top) and t = 2 s (bottom) are shown in 
the x-y plane, which are simulated by hydrodynamic simula¬ 
tions. The size of the snapshot is 256 ^m on each side. The 
color of vortex core represents the directions of mass and spin 
flows = ±l,(Tc, = ±1). Vortices that make annihilation 
pairs have the same symbol shape: red (-I-, -b) and blue (—, —) 
circles, and green (-b, —) and orange ( — ,-b) triangles. 


(<^(^5 O'a') 


V X 'ymass 

(+, +) 

-b 

-b 


-b 

- 

(+; “) 

- 

+ 

(“! +) 

- 

- 


TABLE I. Signs of fz and V x Umass at vortex cores with 

= TIjCTq: = ±1). 


positions of vortex cores are demonstrated in Fig. [2] The 
positions of vortex cores agree with those of maxima of 
V X Uinass- The transverse magnetization and the vortic¬ 
ity of Umass are used to classify vortices into four kinds: 
{cr^,(Ta) = (+,+), (+,-) and (-,+). The com¬ 

bination of (T 0 and Ga is also related to the longitudinal 
magnetization at a vortex core, positive (negative) fz for 
o' 4 , = cfa = —Ca), as mentioned in Sec. Ill Al The sign 
of V X Umass is related to the combination of ctq, and fz 
or simply cr^. Table H] shows the signs of fz and V x Vmass 
at vortex cores for all the combinations of The 

vortices with (-b, -b) and (—, —), which are represented as 
circles in Fig. [2 are a vortex-antivortex pair. Those with 
(-b, —) and (—,-b), which are represented as triangles, are 
another vortex-antivortex pair. 

The number of vortices decreases with time as shown 
in Fig. EKa). In the case of no superfluid flow, which is 
simulated by Eq. m with Wmass = 0 at all times, the 
decay is slower than the other simulations. This sug¬ 
gests that the superfluid flow has an effect to accelerate 
the coarsening dynamics. However, the effect is not very 
simple, which is suggested in Fig. E^b). The dashed line 



FIG. 3. (Golor online) (a) Time dependence of the number 
of vortices simulated by GP equation (GP), hydrodynamic 
equation equation (hydro), and Eq. (1181) in the absence of 
Vmass (no vmass). Each curve is the average of ten simula¬ 
tions. (b) The same data replotted to show how they fit to 
X = (log(l/y) — l)/y (dashed line) by means of the scaling 
Eq. (1^ . 

represents the scaling, Eq. (E51) . The data are fitted to 

t = a[\og{b/Ny) - 1]/Ny + c, ( 36 ) 

where t and are time and the number of vortices, re¬ 
spectively. Note that where L is the system 

size and L = 256 ^m in the simulations. Eor the fitting, 
the data in the range of > 20 are used. The fitting 
parameters a and c correspond to I? jA,A and to, respec¬ 
tively. Parameter 6, which corresponds to is set to 

be a constant value, b = (256/2.4)^. Actually, the core 
size is estimated to be Rc — 2.4 /im in the condensate 
of spin-1 ®^Rb atoms for q/h = 10 Hz. Since the fit¬ 
ting function is modified to be (t — c)b/a = [log(6/A^v) — 
l](6/iVv), the data are plotted as X = [t — c)b/a and 
Y = Ny/b, and they are expected to be on the curve 
X = [log(l/y) — l]/y. The values of fitting parameters 
in Fig. [^b) are (a, c) = (8.7,0.24) in the GP simula¬ 
tion, (8.1,0.050) in the hydrodynamic simulation, and 
(25.4, —0.16) in the absence of Umass- Although the data 
in the absence of t^mass fit to the curve well, those of GP 
and hydrodynamic simulations are very different from the 
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expected scaling. 

It might look strange that the curves of the GP and 
hydrodynamic simulations behave different in Fig. [3] (a), 
although they are similar in Fig.[3(b). Actually, just the 
early-time dynamics is different between GP and hydro- 
dynamic simulations. The given initial states, which are 
unstable, strongly affect the early-time dynamics. After 
the early time, both the GP and hydrodynamic simula¬ 
tions follow a common growth law. Since the growth law 
is not just a power law, the rescaled plots in Fig. |3] (b) 
behave similar even though they look different in the orig¬ 
inal plots. 

The difference of situations between the GP and hydro- 
dynamic simulations and the simulation without Umass is 
twofold. First, the superfluid flow may reduce friction 
(resistivity) in ferromagnetic BECs, which results in the 
faster decay of the number of vortices in the GP and 
hydrodynamic simulations than the simulation without 
■Wmass- Second, there is only one combination of vortex- 
antivortex pair in the absence of superfluid flow, which 
is the same situation in the XY model. In other words, 
when Uinass = 0 , (70 is meaningless, and vortices with the 
opposite signs of aa make a vortex-antivortex pair. By 
contrast, there are two combinations of vortex-antivortex 
pairs [i.e., one is {a^,(7a) = (+,+) and and the 

other is (-I-, — ) and (-I-, —)] in the GP and hydrodynamic 
simulations. 

In order to clarify the reason why the data in the 
above GP and hydrodynamic simulations do not agree 
with the expected scaling, we demonstrate the simula¬ 
tions in special cases where MHVs are limited to two 
kinds that can make a vortex-antivortex pair. MHVs 
of (< 70 , (Tc) = (+,+) are pair-annihilated with those of 
(—,—) but not with the other kinds, (+,—) or (-,-!-). 
If there are only MHVs of (+,+) and (—,—), there is 
only one combination of annihilation pairs, which is the 
same situation as the simulation in the absence of su¬ 
perfluid flow. Then, we can see pure effects of Umass on 
the coarsening dynamics. The simulations with MHVs 
of (+, —) and (—,+) also give the same situation. In 
Fig. mja), the number of vortices decays slightly faster 
than the hydrodynamic simulations, and thus, much 
faster than the simulation in the absence of Umass- This 
fact implies that the coarsening dynamics is accelerated 
by superfluid flow. On the other hand, Fig. HKb) illus¬ 
trates better fitting for the two-kind-vortex data (labeled 
by “pp-|-mm” and “pm-l-mp”) than the four-kind-vortex 
data (labeled by “hydro”). The values of fitting parame¬ 
ters in Fig.|3Kb) are (a, c) = (6.5,-0.030) for “pp-|-mm”, 
and (6.2,-0.030) for “pm-|-mp”. This result indicates 
that the scaling that describes the time dependence of 
vortex density is different from the expected one, Eq (1351) , 
when there are several combinations of vortex-antivortex 
pairs. 



FIG. 4. (Color online) (a) Time dependence of the num¬ 
ber of vortices simulated by hydrodynamic equations in the 
case where there are only two kinds of vortices, (<70, (Ja) = 
(-I-, -b), (—, — ) (labeled with “pp+mm”) or (-)-, —), (—, -I-) (la¬ 
beled with “pm-bmp”). They decay faster than the case where 
four kinds of vortices exist, whose label is “hydro” (the same 
data in Fig. [3] (a)). Each curve is the average of 10 simula¬ 
tions. (b) The same data replotted to show how they fit to 
X = (log(l/y) — 1)/Y (dashed line) by means of the scaling 
Eq. (1^ . 

IV. REVISED GROWTH LAW 

We here consider revising the scaling, and hence the 
growth law, for the case where there are two combina¬ 
tions (groups) of vortex-antivortex pairs. Vortices be¬ 
longing to different groups cannot cause annihilation with 
each other. When the groups of vortex density pi and p 2 
are mixed and coexist in the same space, we rewrite the 
total vortex density ptot = Pi + P 2 as 

Ptot = 2p-2po, (37) 

where p represents a typical vortex density of a group and 
is supposed to obey the original scaling Eq. dSSl), and 2po 
corresponds to the difference between the expected and 
actual vortex densities. Suppose pi = p and p 2 = p—2po, 
where po > 0. This implies that the vortices with density 
pi, which are in the majority, dominate the coarsening 
dynamics. The difference between them pi — p 2 = 2po 













A/fric/A/ 2 , is different between the presence and absence 
of superfluid flow. From Eqs. (iiB, m, and m, 
A/fric/A/2 = —(1 + 3r^)/6 in the presence of superfluid 
flow. When there is no superfluid flow, A/fric/A /2 = 
(1 + r^)A/’{/A/’ 2 , where A// and are given by the 
same winding numbers as the 2D XY model, and thus, 
A/i/A/”^ = —1/2. This means A/fric/A /2 = —(1 + r^)/2 in 
the absence of superfluid flow. Thus, the value of a in 
the simulation in the absence of superfluid flow should 
be about three times larger than that in the GP and hy¬ 
drodynamic simulations. Actually, in Fig. [3](b), a = 8.7 
and 8.1 for the GP and hydrodynamic simulations, re¬ 
spectively, and they are about 1/3 of a = 25.4 for the 
simulation in the absence of superfluid flow. 


FIG. 5. (Color online) The same data of the GP and hydrody¬ 
namic simulations as Fig. replotted to show how they fit to 
X = (log(l/y) — 1)/Y (dashed line) by means of the revised 
scaling Eq. 


is almost independent of time, because the decay rate of 
vortex density is similar to each other if pi ~ p 2 - Even 
if Pi = P 2 in the initial state, difference in vortex density 
often arises in a early time, when the vortex density is 
too high to obey the scaling. The difference is small, 
however, it is the key in the revised scaling. 

Substituting p of Eq. (l37l) into p of Eq. (1^ . we obtain 
a revised equation 


1 log[Pc/(ptot/2-bpo)] - 1 

4 A Ptot/2 -b po 


(38) 


The same data of the GP and hydrodynamic simulations 
as that of Fig. [U which are fitted to Eq. (IMll . are shown 
in Fig. [SJ The data are actually fitted to 


t = a{log[5/(A^v/2 + d)] - l}/(A^v/2 + d) + c (39) 


with b = (256/2.4)^. The fitting parameters a and c cor¬ 
respond to L^/4A and to, respectively, and d corresponds 
to pqL'^. Since the fitting function is modified to be 
{t—c)b/a = {log[6/(Aiv/2-|-(i)] —l}[6/(iVv/2-|-d)], the data 
are plotted as A = {t — c)b/a and Y = {N^/2 + d)/b, and 
they are expected to be on the curve X = [log(l/y) — 
l]/y. The values of fitting parameters in Fig. [5] are 
{a,c,d) = (8.6,-0.059,9.3) and (6.4,—0.19, 6.0) in the 
GP and hydrodynamic simulations, respectively. The 
data are in good agreement with the revised scaling. 


We have considered only MHVs in this paper. The 
coarsening dynamics becomes different and even faster 
in the cases of PCVs and one-component BECs than in 
the case of MHVs. Since some of the assumptions made 
in this paper are invalid for PGVs and one-component 
BECs, the growth laws in those cases should be different 
from that of the XY model or our revised one. We will 
present the study about those cases somewhere else. 

In conclusion, the coarsening dynamics in ferromag¬ 
netic BECs with a positive quadratic Zeeman energy, in 
which magnetic anisotropy is similar to the XY model, 
leads to a different domain growth law from that of the 
XY model. We have proposed a revised growth law es¬ 
pecially for the case where only MHVs exist. When sev¬ 
eral groups of vortex-antivortex pairs coexist in the same 
space, the difference in vortex densities of them leads 
to the revised growth law. In the absence of the super- 
fluid flow, where there is only one combination of vortex- 
antivortex pairs, the growth law is the same as that of 
the XY model, and the coarsening dynamics is slower 
than in the presence of the flow. The effect of the su¬ 
perfluid flow is not only accelerating domain growth but 
also producing several combinations of vortex-antivortex 
pairs. 
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